#How Do Observers Assess Resolve?
#July 7, 2018
#BJPS Replication 4.R

#This file contains the replication code needed to replicate the analysis from Supplementary Appendix §1 - §1.1 (model diagnostics, randomization checks, dropping unusual dyadic combinations)

#To generate the necessary dataframes, either run BJPS Replication 1.R, or load the following saved objects:
library(here) #Avoids needing to setwd()
dat3 <- get(load(file="Resolve Conjoint Cleaned No US 070718.RData"))
dat2 <- get(load(file="Resolve Conjoint Cleaned 070718.RData"))
dat_id <- get(load(file="Resolve Demographics 070718.RData"))

#Load parallel processing and recode functions in case you haven't run BJPS Replication 2.R - if not, you can skip to line 58
library(snow)
cl <- makeCluster(8,"SOCK")
clusterBootS2 <- function(data){
	i <- sample(unique(data$ID_numeric),length(unique(data$ID_numeric)),replace=TRUE)
	row.nums <- NULL
	for (j in 1:length(i)){
		row.nums <- c(row.nums, which(data$ID_numeric==i[j]))
	}
	return(data[row.nums,])
}

#Functions to recode our variables
recode2 <- function(variable, reverse=FALSE, maxVal, minVal, binarize=FALSE, x=NA, x2){
	if (is.factor(variable)){variable <- as.numeric(variable)}
	if (is.character(variable)){variable <- as.numeric(variable)}	
	if (missing(maxVal)){maxVal <- max(variable, na.rm=TRUE)}
	if (missing(minVal)){minVal <- min(variable, na.rm=TRUE)}
	variable[variable > maxVal] <- NA
	variable[variable < minVal] <- NA
	if (reverse){
		temp <- variable
		for (j in 1:maxVal){
			temp[which(variable==j)] <- maxVal-j+1
		}
		return(as.numeric(temp))
	}
	if (binarize){
		if (is.na(x)){stop("Specify the value to dichotomize on")}
		temp <- variable
		if (missing(x2)){
			i <- which(variable == x)
			k <- which(variable < x | variable > x)
		}
		else{
			i <- which(variable >= x & variable <= x2)
			k <- which(variable < x | variable > x2)
		}
		temp[i] <- 1
		temp[k] <- 0
		return(temp)
	}
	return(as.numeric(variable))
}

##### Testing stability and carryover effects assumption: how do AMCEs change over time?

#Reference model
mod.conj2 <- lm(Chosen ~ RegimeType + Capabilities + Stakes + NewLeader + milService + maleLeader + actor + prevActLdr + prevOpponent + prevRole + curBehavior, data=dat3)

stability.results <- matrix(NA, nrow=length(coef(mod.conj2)), ncol=8)
rownames(stability.results) <- names(coef(mod.conj2))
for (i in 1:8){
	temp <- lm(Chosen ~ RegimeType + Capabilities + Stakes + NewLeader + milService + maleLeader + actor + prevActLdr + prevOpponent + prevRole + curBehavior, data=subset(dat3, dat3$Round==i))
	stability.results[,i] <- coef(temp)
}

round(stability.results, digits=3) #Results differ in sixth choice task

#Estimate results dropping round 6
set.seed(43215)
dat_6 <- subset(dat3, dat3$Round!=6)
bootConjoint_6 <- function(...){
	temp <- clusterBootS2(dat_6)
	mod.temp <- lm(Chosen ~ Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("dat_6", "bootConjoint_6", "clusterBootS2"))

system.time(boot.6 <- parSapply(cl, rep(1,1500), bootConjoint_6)) 

#Extract quantities of interest
plot.mat6 <- cbind(apply(boot.6, 1, mean), apply(boot.6, 1, quantile, c(0.025)), apply(boot.6, 1, quantile, c(0.05)), apply(boot.6, 1, quantile, c(0.95)), apply(boot.6,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.mat6b <- rbind(rep(0,5), plot.mat6[1,], rep(0,5), plot.mat6[2,], rep(0,5), plot.mat6[3:4,], rep(0,5), plot.mat6[5,], rep(0,5), plot.mat6[6,], rep(0,5), plot.mat6[c(8,7),], rep(0,5), plot.mat6[9,], rep(0,5), plot.mat6[10,], rep(0,5), plot.mat6[11,], rep(0,5), plot.mat6[c(14,13,12),],rep(0,5), plot.mat6[15:16,])

textLab <- c("Low Capabilities", "High Capabilities", "Low Stakes", "High Stakes", "Dictatorship", "Democracy", "Mixed",  "Ally", "Adversary", "Established Leader", "New Leader", "No Military Service", "Some Military Service", "Long Military Service", "Female Leader", "Male Leader", "Target", "Initiator", "Against ally", "Against adversary", "Different leader, stood firm", "Same leader, stood firm", "Same leader, backed down", "Different leader, backed down", "Nothing", "Mobilized troops", "Public threat")

##### Figure 1: Dropping the sixth choice task #####
# Clean up labels in post-production
#####
dev.new(height=9, width=6.5)
par(mar=c(5.1,1.1,4.1,0), oma=c(0,4,0,0))
plot(plot.mat2b[,1], nrow(plot.mat2b):1, type="n", axes=FALSE, xlab="Average marginal component effect (AMCE)", ylab="",xlim=c(-0.3,0.2))

for (i in 1:nrow(plot.mat2b)){
	points(plot.mat2b[i,1], nrow(plot.mat2b)-(i-1), pch=16)
	lines(plot.mat2b[i,c(2,5)], rep(nrow(plot.mat2b)-(i-1),2))
	lines(plot.mat2b[i,c(3,4)], rep(nrow(plot.mat2b)-(i-1),2), lwd=2)		
	text(-0.3, nrow(plot.mat2b)-(i-1), textLab[i], cex=0.7)			
	}
abline(v=0)
axis(1)

##### Testing profile order assumption: do results vary when attribute appears for country A vs country B?

#Bootstrap separate models for Country A vs Country B

datA <- subset(dat3, dat3$Country=="A")
datB <- subset(dat3, dat3$Country=="B")

bootConjointA <- function(...){
	temp <- clusterBootS2(datA)
	mod.temp <- lm(Chosen ~ Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=temp)
	return(coef(mod.temp))
}

bootConjointB <- function(...){
	temp <- clusterBootS2(datB)
	mod.temp <- lm(Chosen ~ Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("datA", "datB", "bootConjointA", "bootConjointB","clusterBootS2"))

system.time(boot.fullA <- parSapply(cl, rep(1,1500), bootConjointA)) #57 seconds!
system.time(boot.fullB <- parSapply(cl, rep(1,1500), bootConjointB)) #57 seconds!

#Extract quantities of interest
plot.fullA <- cbind(apply(boot.fullA, 1, mean), apply(boot.fullA, 1, quantile, c(0.025)), apply(boot.fullA, 1, quantile, c(0.05)), apply(boot.fullA, 1, quantile, c(0.95)), apply(boot.fullA,1,quantile, c(0.975)))[-1,]
plot.fullB <- cbind(apply(boot.fullB, 1, mean), apply(boot.fullB, 1, quantile, c(0.025)), apply(boot.fullB, 1, quantile, c(0.05)), apply(boot.fullB, 1, quantile, c(0.95)), apply(boot.fullB,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.fullAb <- rbind(rep(0,5), plot.fullA[1,], rep(0,5), plot.fullA[2,], rep(0,5), plot.fullA[3:4,], rep(0,5), plot.fullA[5,], rep(0,5), plot.fullA[6,], rep(0,5), plot.fullA[c(8,7),], rep(0,5), plot.fullA[9,], rep(0,5), plot.fullA[10,], rep(0,5), plot.fullA[11,], rep(0,5), plot.fullA[c(14,13,12),],rep(0,5), plot.fullA[15:16,])
plot.fullBb <- rbind(rep(0,5), plot.fullB[1,], rep(0,5), plot.fullB[2,], rep(0,5), plot.fullB[3:4,], rep(0,5), plot.fullB[5,], rep(0,5), plot.fullB[6,], rep(0,5), plot.fullB[c(8,7),], rep(0,5), plot.fullB[9,], rep(0,5), plot.fullB[10,], rep(0,5), plot.fullB[11,], rep(0,5), plot.fullB[c(14,13,12),],rep(0,5), plot.fullB[15:16,])

#Labels
textLab <- c("Low Capabilities", "High Capabilities", "Low Stakes", "High Stakes", "Dictatorship", "Democracy", "Mixed",  "Ally", "Adversary", "Established Leader", "New Leader", "No Military Service", "Some Military Service", "Long Military Service", "Female Leader", "Male Leader", "Target", "Initiator", "Against ally", "Against adversary", "Different leader, stood firm", "Same leader, stood firm", "Same leader, backed down", "Different leader, backed down", "Nothing", "Mobilized troops", "Public threat")

##### Figure 2: testing the profile order assumption ####
# Fix labels in post-production

dev.new(height=9, width=6.5)
par(mar=c(5.1,1.1,4.1,0), oma=c(0,4,0,0))
plot(plot.fullAb[,1], nrow(plot.fullAb):1, type="n", axes=FALSE, xlab="Average marginal component effect (AMCE)", ylab="",xlim=c(-0.3,0.2))

for (i in 1:nrow(plot.fullAb)){
	points(plot.fullAb[i,1], nrow(plot.fullAb)-(i-1), pch=16)
	points(plot.fullBb[i,1], nrow(plot.fullBb)-(i-1)-0.25, pch=16, col="grey")
	lines(plot.fullAb[i,c(2,5)], rep(nrow(plot.fullAb)-(i-1),2))
	lines(plot.fullBb[i,c(2,5)], rep(nrow(plot.fullBb)-(i-1)-0.25,2), col="grey")
	lines(plot.fullAb[i,c(3,4)], rep(nrow(plot.fullAb)-(i-1),2), lwd=2)	
	lines(plot.fullBb[i,c(3,4)], rep(nrow(plot.fullAb)-(i-1)-0.25,2), lwd=2, col="grey")		
	text(-0.3, nrow(plot.fullAb)-(i-1), textLab[i], cex=0.7)			
	}
abline(v=0)
axis(1)

###### Randomization check

#First, merge demographics into main dataset
dat.merge <- merge(dat3, dat_id, by="ID_numeric")

#Now estimate a series of logits (for factors with 2 levels) or multinomial logits (for factors with >2 levels) where each treatment assignment is modeled as a function of a set of demographic variables


bootRandomCheck <- function(...){
	temp <- clusterBootS2(dat.merge)
	mod.temp <- glm(as.factor(Capabilities) ~ male + age + educ + partyID + mi1, data=temp, family=binomial)
	return(coef(mod.temp))
}

clusterExport(cl, list("dat.merge", "bootRandomCheck"))

system.time(r1 <- parSapply(cl, rep(1,1000), bootRandomCheck)) 

r1.qi <- round(cbind(apply(r1, 1, mean), apply(r1, 1, quantile, c(0.025)), apply(r1,1,quantile, c(0.975))), digits=3) #None significant at p < 0.05

bootRandomCheck <- function(...){
	temp <- clusterBootS2(dat.merge)
	mod.temp <- glm(as.factor(Stakes) ~ male + age + educ + partyID + mi1, data=temp, family=binomial)
	return(coef(mod.temp))
}

clusterExport(cl, list("bootRandomCheck"))

system.time(r2 <- parSapply(cl, rep(1,1000), bootRandomCheck))

r2.qi <- round(cbind(apply(r2, 1, mean), apply(r2, 1, quantile, c(0.025)), apply(r2,1,quantile, c(0.975))), digits=3) 

bootRandomCheck <- function(...){
	temp <- clusterBootS2(dat.merge)
	mod.temp <- glm(as.factor(NewLeader) ~ male + age + educ + partyID + mi1, data=temp, family=binomial)
	return(coef(mod.temp))
}

clusterExport(cl, list("bootRandomCheck"))
system.time(r3 <- parSapply(cl, rep(1,1000), bootRandomCheck)) 

r3.qi <- round(cbind(apply(r3, 1, mean), apply(r3, 1, quantile, c(0.025)), apply(r3,1,quantile, c(0.975))), digits=3) 

bootRandomCheck <- function(...){
	temp <- clusterBootS2(dat.merge)
	mod.temp <- glm(as.factor(maleLeader) ~ male + age + educ + partyID + mi1, data=temp, family=binomial)
	return(coef(mod.temp))
}

clusterExport(cl, list("bootRandomCheck"))
system.time(r4 <- parSapply(cl, rep(1,1000), bootRandomCheck)) 

r4.qi <- round(cbind(apply(r4, 1, mean), apply(r4, 1, quantile, c(0.025)), apply(r4,1,quantile, c(0.975))), digits=3) 

bootRandomCheck <- function(...){
	temp <- clusterBootS2(dat.merge)
	mod.temp <- glm(as.factor(actor) ~ male + age + educ + partyID + mi1, data=temp, family=binomial)
	return(coef(mod.temp))
}

clusterExport(cl, list("bootRandomCheck"))

system.time(r5 <- parSapply(cl, rep(1,1000), bootRandomCheck)) 

r5.qi <- round(cbind(apply(r5, 1, mean), apply(r5, 1, quantile, c(0.025)), apply(r5,1,quantile, c(0.975))), digits=3) #Partisanship significant?

bootRandomCheck <- function(...){
	temp <- clusterBootS2(dat.merge)
	mod.temp <- glm(as.factor(prevOpponent) ~ male + age + educ + partyID + mi1, data=temp, family=binomial)
	return(coef(mod.temp))
}

clusterExport(cl, list("bootRandomCheck"))

system.time(r6 <- parSapply(cl, rep(1,1000), bootRandomCheck)) 

r6.qi <- round(cbind(apply(r6, 1, mean), apply(r6, 1, quantile, c(0.025)), apply(r6,1,quantile, c(0.975))), digits=3) #None significant at p < 0.05

bootRandomCheck <- function(...){
	temp <- clusterBootS2(dat.merge)
	mod.temp <- glm(as.factor(prevRole) ~ male + age + educ + partyID + mi1, data=temp, family=binomial)
	return(coef(mod.temp))
}

clusterExport(cl, list("bootRandomCheck"))

system.time(r7 <- parSapply(cl, rep(1,1000), bootRandomCheck))

r7.qi <- round(cbind(apply(r7, 1, mean), apply(r7, 1, quantile, c(0.025)), apply(r7,1,quantile, c(0.975))), digits=3) 

bootRandomCheck <- function(...){
	temp <- clusterBootS2(dat.merge)
	mod.temp <- glm(as.factor(prevAct) ~ male + age + educ + partyID + mi1, data=temp, family=binomial)
	return(coef(mod.temp))
}

clusterExport(cl, list("bootRandomCheck"))

system.time(r8 <- parSapply(cl, rep(1,1000), bootRandomCheck)) 

r8.qi <- round(cbind(apply(r8, 1, mean), apply(r8, 1, quantile, c(0.025)), apply(r8,1,quantile, c(0.975))), digits=3) 

bootRandomCheck <- function(...){
	temp <- clusterBootS2(dat.merge)
	mod.temp <- glm(as.factor(prevLdr) ~ male + age + educ + partyID + mi1, data=temp, family=binomial)
	return(coef(mod.temp))
}

clusterExport(cl, list("bootRandomCheck"))

system.time(r9 <- parSapply(cl, rep(1,1000), bootRandomCheck)) 

r9.qi <- round(cbind(apply(r9, 1, mean), apply(r9, 1, quantile, c(0.025)), apply(r9,1,quantile, c(0.975))), digits=3)

#Need multinomial for regime type, mil service, current behavior

library(nnet)
bootRandomCheck <- function(...){
	temp <- clusterBootS2(dat.merge)
	mod.temp <- multinom(as.factor(RegimeType) ~ male + age + educ + partyID + mi1, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("bootRandomCheck", "multinom"))

system.time(r10 <- parSapply(cl, rep(1,1000), bootRandomCheck))

r10.qi <- round(cbind(apply(r10, 1, mean), apply(r10, 1, quantile, c(0.025)), apply(r10,1,quantile, c(0.975))), digits=3) 

bootRandomCheck <- function(...){
	temp <- clusterBootS2(dat.merge)
	mod.temp <- multinom(as.factor(milService) ~ male + age + educ + partyID + mi1, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("bootRandomCheck", "multinom"))

system.time(r11 <- parSapply(cl, rep(1,1000), bootRandomCheck)) #101 seconds

r11.qi <- round(cbind(apply(r11, 1, mean), apply(r11, 1, quantile, c(0.025)), apply(r11,1,quantile, c(0.975))), digits=3) 

bootRandomCheck <- function(...){
	temp <- clusterBootS2(dat.merge)
	mod.temp <- multinom(as.factor(curBehavior) ~ male + age + educ + partyID + mi1, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("bootRandomCheck", "multinom"))

system.time(r12 <- parSapply(cl, rep(1,1000), bootRandomCheck)) 

r12.qi <- round(cbind(apply(r12, 1, mean), apply(r12, 1, quantile, c(0.025)), apply(r12,1,quantile, c(0.975))), digits=3) 

##### Table 1: Randomization Check #####

#To generate table with results, only do r1-9 first, and then separate the MNLs by DV (slightly jerry-rigged and not as elegant as it should be, don't judge)

library(stargazer)
tabRCheck <- cbind(r1.qi[,2:3], r2.qi[,2:3], r3.qi[,2:3],r4.qi[,2:3], r5.qi[,2:3], r6.qi[,2:3], r7.qi[,2:3], r8.qi[,2:3], r9.qi[,2:3], r10.qi[seq(1,12,by=2),2:3], r10.qi[seq(2,12,by=2),2:3], r11.qi[seq(1,12,by=2),2:3], r11.qi[seq(2,12,by=2),2:3], r12.qi[seq(1,12,by=2),2:3], r12.qi[seq(2,12,by=2),2:3])

tabRCheck2 <- matrix(NA, nrow=nrow(tabRCheck), ncol=ncol(tabRCheck)/2)
rownames(tabRCheck2) <- c("(Intercept)", "Male", "Age", "Education", "Party ID", "Mil Assert")
colnames(tabRCheck2) <- c("High capabilities", "High stakes", "New leader", "Male leader", "Adversary", "Against adversary", "Initiator", "Backed down", "Same leader", "Democracy", "Mixed", "High service", "Some service", "Mobilized troops", "Public threat")
for (i in 1:ncol(tabRCheck2)){
	for (j in 1:nrow(tabRCheck2)){
		tabRCheck2[j,i] <- paste("[",tabRCheck[j,(2*i)-1], ", ", tabRCheck[j,(2*i)], "]", sep="")
	}
}

#Assemble table below
stargazer(tabRCheck2[,1:4])
stargazer(tabRCheck2[,5:8])
stargazer(tabRCheck2[,13:15])

##### Attribute order tests (e.g. whether it matters whether an attribute was in the first row or the seventh)

#To retrieve order information, return to the original raw dataset
data <- get(load("Resolve Conjoint Raw 070718.RData"))

#Create miniature version of dataset with ordering information; note that the ordering will be the same across all tasks per respondent, but vary between respondents
datOrder <- data.frame(ID_numeric=unique(dat3$ID_numeric), ID=unique(dat3$ID), RegimeTypeO=NA, CapabilitiesO=NA, StakesO=NA, LeaderO=NA, actorO=NA, prevBehaviorO=NA, curBehaviorO=NA)

for (i in 1:nrow(datOrder)){
	loc <- which(data$V1 ==datOrder$ID[i])
	if(!is.na(loc)){
	ordering <- data[loc, seq(11,23,by=2)] #Just extracts ordering from country A, since it should be the same for both
	datOrder$RegimeTypeO[i] <- which(ordering=="Government type")
	datOrder$CapabilitiesO[i] <- which(ordering=="Military capabilities")
	datOrder$StakesO[i] <- which(ordering=="Interests in the dispute")
	datOrder$LeaderO[i] <-  which(ordering=="Leader background")
	datOrder$actorO[i] <- which(ordering=="Foreign relations")
	datOrder$prevBehaviorO[i] <- which(ordering=="Previous behavior in international disputes")
	datOrder$curBehaviorO[i] <- which(ordering=="Current behavior")
	}
}

#Merge in order variables to main dataset
dat.merge2 <- merge(datOrder, dat3, by="ID_numeric")

#Now, test whether the order of attributes makes a difference

dat.merge2$OCap1 <- recode2(dat.merge2$CapabilitiesO, binarize=TRUE, x=1)
dat.merge2$OCap2 <- recode2(dat.merge2$CapabilitiesO, binarize=TRUE, x=2)
dat.merge2$OCap3 <- recode2(dat.merge2$CapabilitiesO, binarize=TRUE, x=3)
dat.merge2$OCap4 <- recode2(dat.merge2$CapabilitiesO, binarize=TRUE, x=4)
dat.merge2$OCap5 <- recode2(dat.merge2$CapabilitiesO, binarize=TRUE, x=5)
dat.merge2$OCap6 <- recode2(dat.merge2$CapabilitiesO, binarize=TRUE, x=6)
dat.merge2$OCap7 <- recode2(dat.merge2$CapabilitiesO, binarize=TRUE, x=7)

bootOrderCheck <- function(...){
	temp <- clusterBootS2(dat.merge2)
	mod.temp <- lm(Chosen ~ Capabilities * (OCap2 + OCap3 + OCap4 + OCap5 + OCap6 + OCap7), data=temp)
	x <- coef(mod.temp)[c(2,9,10,11,12,13,14)] #Extract coefficient for treatment of interest, plus its interaction terms
	x[2:7] <- x[2:7]+x[1] #Add treatment coefficient to each interaction coefficient to get the conditional effect of the treatment at each attribute ordering
	return(x)
}

clusterExport(cl, list("dat.merge2", "bootOrderCheck", "clusterBootS2"))

system.time(o1 <- parSapply(cl, rep(1,1000), bootOrderCheck)) 

dat.merge2$OStakes1 <- recode2(dat.merge2$StakesO, binarize=TRUE, x=1)
dat.merge2$OStakes2 <- recode2(dat.merge2$StakesO, binarize=TRUE, x=2)
dat.merge2$OStakes3 <- recode2(dat.merge2$StakesO, binarize=TRUE, x=3)
dat.merge2$OStakes4 <- recode2(dat.merge2$StakesO, binarize=TRUE, x=4)
dat.merge2$OStakes5 <- recode2(dat.merge2$StakesO, binarize=TRUE, x=5)
dat.merge2$OStakes6 <- recode2(dat.merge2$StakesO, binarize=TRUE, x=6)
dat.merge2$OStakes7 <- recode2(dat.merge2$StakesO, binarize=TRUE, x=7)

bootOrderCheck <- function(...){
	temp <- clusterBootS2(dat.merge2)
	mod.temp <- lm(Chosen ~ Stakes * (OStakes2 + OStakes3 + OStakes4 + OStakes5 + OStakes6 + OStakes7), data=temp)
	x <- coef(mod.temp)[c(2,9,10,11,12,13,14)]
	x[2:7] <- x[2:7]+x[1]
	return(x)
}

clusterExport(cl, list("dat.merge2", "bootOrderCheck", "clusterBootS2"))

system.time(o2 <- parSapply(cl, rep(1,1000), bootOrderCheck))

dat.merge2$ORegimeType1 <- recode2(dat.merge2$RegimeTypeO, binarize=TRUE, x=1)
dat.merge2$ORegimeType2 <- recode2(dat.merge2$RegimeTypeO, binarize=TRUE, x=2)
dat.merge2$ORegimeType3 <- recode2(dat.merge2$RegimeTypeO, binarize=TRUE, x=3)
dat.merge2$ORegimeType4 <- recode2(dat.merge2$RegimeTypeO, binarize=TRUE, x=4)
dat.merge2$ORegimeType5 <- recode2(dat.merge2$RegimeTypeO, binarize=TRUE, x=5)
dat.merge2$ORegimeType6 <- recode2(dat.merge2$RegimeTypeO, binarize=TRUE, x=6)
dat.merge2$ORegimeType7 <- recode2(dat.merge2$RegimeTypeO, binarize=TRUE, x=7)

bootOrderCheck <- function(...){
	temp <- clusterBootS2(dat.merge2)
	mod.temp <- lm(Chosen ~ RegimeType * (ORegimeType2 + ORegimeType3 + ORegimeType4 + ORegimeType5 + ORegimeType6 + ORegimeType7), data=temp)
	x <- coef(mod.temp)[c(2,10,12,14,16,18,20,3,11,13,15,17,19,21)]
	x[2:7] <- x[2:7]+x[1]
	x[9:14] <- x[9:14] + x[8]
	return(x)
}

clusterExport(cl, list("dat.merge2", "bootOrderCheck", "clusterBootS2"))

system.time(o3 <- parSapply(cl, rep(1,1000), bootOrderCheck)) 

dat.merge2$OLeader1 <- recode2(dat.merge2$LeaderO, binarize=TRUE, x=1)
dat.merge2$OLeader2 <- recode2(dat.merge2$LeaderO, binarize=TRUE, x=2)
dat.merge2$OLeader3 <- recode2(dat.merge2$LeaderO, binarize=TRUE, x=3)
dat.merge2$OLeader4 <- recode2(dat.merge2$LeaderO, binarize=TRUE, x=4)
dat.merge2$OLeader5 <- recode2(dat.merge2$LeaderO, binarize=TRUE, x=5)
dat.merge2$OLeader6 <- recode2(dat.merge2$LeaderO, binarize=TRUE, x=6)
dat.merge2$OLeader7 <- recode2(dat.merge2$LeaderO, binarize=TRUE, x=7)

bootOrderCheck <- function(...){
	temp <- clusterBootS2(dat.merge2)
	mod.temp <- lm(Chosen ~ NewLeader * (OLeader2 + OLeader3 + OLeader4 + OLeader5 + OLeader6 + OLeader7), data=temp)
	x <- coef(mod.temp)[c(2,9,10,11,12,13,14)]
	x[2:7] <- x[2:7]+x[1]
	return(x)
}

clusterExport(cl, list("dat.merge2", "bootOrderCheck", "clusterBootS2"))

system.time(o4 <- parSapply(cl, rep(1,1000), bootOrderCheck))

bootOrderCheck <- function(...){
	temp <- clusterBootS2(dat.merge2)
	mod.temp <- lm(Chosen ~ milService * (OLeader2 + OLeader3 + OLeader4 + OLeader5 + OLeader6 + OLeader7), data=temp)
	x <- coef(mod.temp)[c(2,10,12,14,16,18,20,3,11,13,15,17,19,21)]
	x[2:7] <- x[2:7]+x[1]
	x[9:14] <- x[9:14] + x[8]
	return(x)
}

clusterExport(cl, list("dat.merge2", "bootOrderCheck", "clusterBootS2"))

system.time(o5 <- parSapply(cl, rep(1,1000), bootOrderCheck)) 

bootOrderCheck <- function(...){
	temp <- clusterBootS2(dat.merge2)
	mod.temp <- lm(Chosen ~ maleLeader * (OLeader2 + OLeader3 + OLeader4 + OLeader5 + OLeader6 + OLeader7), data=temp)
	x <- coef(mod.temp)[c(2,9,10,11,12,13,14)]
	x[2:7] <- x[2:7]+x[1]
	return(x)
}

clusterExport(cl, list("dat.merge2", "bootOrderCheck", "clusterBootS2"))

system.time(o6 <- parSapply(cl, rep(1,1000), bootOrderCheck)) 

dat.merge2$OActor1 <- recode2(dat.merge2$actorO, binarize=TRUE, x=1)
dat.merge2$OActor2 <- recode2(dat.merge2$actorO, binarize=TRUE, x=2)
dat.merge2$OActor3 <- recode2(dat.merge2$actorO, binarize=TRUE, x=3)
dat.merge2$OActor4 <- recode2(dat.merge2$actorO, binarize=TRUE, x=4)
dat.merge2$OActor5 <- recode2(dat.merge2$actorO, binarize=TRUE, x=5)
dat.merge2$OActor6 <- recode2(dat.merge2$actorO, binarize=TRUE, x=6)
dat.merge2$OActor7 <- recode2(dat.merge2$actorO, binarize=TRUE, x=7)

bootOrderCheck <- function(...){
	temp <- clusterBootS2(dat.merge2)
	mod.temp <- lm(Chosen ~ actor * (OActor2 + OActor3 + OActor4 + OActor5 + OActor6 + OActor7), data=temp)
	x <- coef(mod.temp)[c(2,9,10,11,12,13,14)]
	x[2:7] <- x[2:7]+x[1]
	return(x)
}

clusterExport(cl, list("dat.merge2", "bootOrderCheck", "clusterBootS2"))

system.time(o7 <- parSapply(cl, rep(1,1000), bootOrderCheck)) 

dat.merge2$OprevBehavior1 <- recode2(dat.merge2$prevBehaviorO, binarize=TRUE, x=1)
dat.merge2$OprevBehavior2 <- recode2(dat.merge2$prevBehaviorO, binarize=TRUE, x=2)
dat.merge2$OprevBehavior3 <- recode2(dat.merge2$prevBehaviorO, binarize=TRUE, x=3)
dat.merge2$OprevBehavior4 <- recode2(dat.merge2$prevBehaviorO, binarize=TRUE, x=4)
dat.merge2$OprevBehavior5 <- recode2(dat.merge2$prevBehaviorO, binarize=TRUE, x=5)
dat.merge2$OprevBehavior6 <- recode2(dat.merge2$prevBehaviorO, binarize=TRUE, x=6)
dat.merge2$OprevBehavior7 <- recode2(dat.merge2$prevBehaviorO, binarize=TRUE, x=7)

bootOrderCheck <- function(...){
	temp <- clusterBootS2(dat.merge2)
	mod.temp <- lm(Chosen ~ prevOpponent * (OprevBehavior2 + OprevBehavior3 + OprevBehavior4 + OprevBehavior5 + OprevBehavior6 + OprevBehavior7), data=temp)
	x <- coef(mod.temp)[c(2,9,10,11,12,13,14)]
	x[2:7] <- x[2:7]+x[1]
	return(x)
}

clusterExport(cl, list("dat.merge2", "bootOrderCheck", "clusterBootS2"))

system.time(o8 <- parSapply(cl, rep(1,1000), bootOrderCheck)) 


bootOrderCheck <- function(...){
	temp <- clusterBootS2(dat.merge2)
	mod.temp <- lm(Chosen ~ prevRole * (OprevBehavior2 + OprevBehavior3 + OprevBehavior4 + OprevBehavior5 + OprevBehavior6 + OprevBehavior7), data=temp)
	x <- coef(mod.temp)[c(2,9,10,11,12,13,14)]
	x[2:7] <- x[2:7]+x[1]
	return(x)
}

clusterExport(cl, list("dat.merge2", "bootOrderCheck", "clusterBootS2"))

system.time(o9 <- parSapply(cl, rep(1,1000), bootOrderCheck)) 

bootOrderCheck <- function(...){
	temp <- clusterBootS2(dat.merge2)
	mod.temp <- lm(Chosen ~ prevAct * (OprevBehavior2 + OprevBehavior3 + OprevBehavior4 + OprevBehavior5 + OprevBehavior6 + OprevBehavior7), data=temp)
	x <- coef(mod.temp)[c(2,9,10,11,12,13,14)]
	x[2:7] <- x[2:7]+x[1]
	return(x)
}

clusterExport(cl, list("dat.merge2", "bootOrderCheck", "clusterBootS2"))

system.time(o10 <- parSapply(cl, rep(1,1000), bootOrderCheck)) 

bootOrderCheck <- function(...){
	temp <- clusterBootS2(dat.merge2)
	mod.temp <- lm(Chosen ~ prevLdr * (OprevBehavior2 + OprevBehavior3 + OprevBehavior4 + OprevBehavior5 + OprevBehavior6 + OprevBehavior7), data=temp)
	x <- coef(mod.temp)[c(2,9,10,11,12,13,14)]
	x[2:7] <- x[2:7]+x[1]
	return(x)
}

clusterExport(cl, list("dat.merge2", "bootOrderCheck", "clusterBootS2"))

system.time(o11 <- parSapply(cl, rep(1,1000), bootOrderCheck)) 

dat.merge2$OcurBehavior1 <- recode2(dat.merge2$curBehaviorO, binarize=TRUE, x=1)
dat.merge2$OcurBehavior2 <- recode2(dat.merge2$curBehaviorO, binarize=TRUE, x=2)
dat.merge2$OcurBehavior3 <- recode2(dat.merge2$curBehaviorO, binarize=TRUE, x=3)
dat.merge2$OcurBehavior4 <- recode2(dat.merge2$curBehaviorO, binarize=TRUE, x=4)
dat.merge2$OcurBehavior5 <- recode2(dat.merge2$curBehaviorO, binarize=TRUE, x=5)
dat.merge2$OcurBehavior6 <- recode2(dat.merge2$curBehaviorO, binarize=TRUE, x=6)
dat.merge2$OcurBehavior7 <- recode2(dat.merge2$curBehaviorO, binarize=TRUE, x=7)

bootOrderCheck <- function(...){
	temp <- clusterBootS2(dat.merge2)
	mod.temp <- lm(Chosen ~ curBehavior * (OcurBehavior2 + OcurBehavior3 + OcurBehavior4 + OcurBehavior5 + OcurBehavior6 + OcurBehavior7), data=temp)
	x <- coef(mod.temp)[c(2,10,12,14,16,18,20,3,11,13,15,17,19,21)]
	x[2:7] <- x[2:7]+x[1]
	x[9:14] <- x[9:14] + x[8]
	return(x)

}

clusterExport(cl, list("dat.merge2", "bootOrderCheck", "clusterBootS2"))

system.time(o12 <- parSapply(cl, rep(1,1000), bootOrderCheck)) 

#Extract quantities of interest
o1.r <- cbind(apply(o1, 1, mean), apply(o1, 1, quantile, c(0.025)), apply(o1,1,quantile, c(0.975)))
o2.r <- cbind(apply(o2, 1, mean), apply(o2, 1, quantile, c(0.025)), apply(o2,1,quantile, c(0.975)))
o3.r <- cbind(apply(o3, 1, mean), apply(o3, 1, quantile, c(0.025)), apply(o3,1,quantile, c(0.975)))
o4.r <- cbind(apply(o4, 1, mean), apply(o4, 1, quantile, c(0.025)), apply(o4,1,quantile, c(0.975)))
o5.r <- cbind(apply(o5, 1, mean), apply(o5, 1, quantile, c(0.025)), apply(o5,1,quantile, c(0.975)))
o6.r <- cbind(apply(o6, 1, mean), apply(o6, 1, quantile, c(0.025)), apply(o6,1,quantile, c(0.975)))
o7.r <- cbind(apply(o7, 1, mean), apply(o7, 1, quantile, c(0.025)), apply(o7,1,quantile, c(0.975)))
o8.r <- cbind(apply(o8, 1, mean), apply(o8, 1, quantile, c(0.025)), apply(o8,1,quantile, c(0.975)))
o9.r <- cbind(apply(o9, 1, mean), apply(o9, 1, quantile, c(0.025)), apply(o9,1,quantile, c(0.975)))
o10.r <- cbind(apply(o10, 1, mean), apply(o10, 1, quantile, c(0.025)), apply(o10,1,quantile, c(0.975)))
o11.r <- cbind(apply(o11, 1, mean), apply(o11, 1, quantile, c(0.025)), apply(o11,1,quantile, c(0.975)))
o12.r <- cbind(apply(o12, 1, mean), apply(o12, 1, quantile, c(0.025)), apply(o12,1,quantile, c(0.975)))

textLab <- c("AMCE 1", "AMCE 2", "AMCE 3", "AMCE 4", "AMCE 5", "AMCE 6", "AMCE 7")

##### Figure 3: do treatment effects systematically vary by attribute order?
dev.new(height=9,width=12)
par(mar=c(5.1,1.1,4.1,0), oma=c(0,0,0,0), mfrow=c(3,4))
plot(o1.r[,1], 1:7, type="n", xlab="AMCE", ylab="", axes=FALSE, xlim=c(-0.35,0.3), ylim=c(0.5,7), main= "Capabilities")
for (i in 1:nrow(o1.r)){
	points(o1.r[i,1], nrow(o1.r)-(i-1), pch=16)
	lines(o1.r[i,c(2,3)], rep(nrow(o1.r)-(i-1),2))
	text(-0.3, nrow(o1.r)-(i-1), textLab[i], cex=0.8)
}
abline(v=0)
axis(1)

plot(o2.r[,1], 1:7, type="n", xlab="AMCE", ylab="", axes=FALSE, xlim=c(-0.35,0.3), ylim=c(0.5,7), main= "Stakes")
for (i in 1:nrow(o2.r)){
	points(o2.r[i,1], nrow(o2.r)-(i-1), pch=16)
	lines(o2.r[i,c(2,3)], rep(nrow(o2.r)-(i-1),2))
	text(-0.3, nrow(o2.r)-(i-1), textLab[i], cex=0.8)
}
abline(v=0)
axis(1)

plot(o3.r[,1], 1:14, type="n", xlab="AMCE", ylab="", axes=FALSE, xlim=c(-0.35,0.3), ylim=c(0.5,14), main= "Regime type")
for (i in 1:nrow(o3.r)){
	points(o3.r[i,1], nrow(o3.r)-(i-1), pch=16)
	lines(o3.r[i,c(2,3)], rep(nrow(o3.r)-(i-1),2))
	text(-0.3, nrow(o3.r)-(i-1), textLab[ifelse(i<8,i,i - 7)], cex=0.8)
}
abline(v=0)
axis(1)

plot(o4.r[,1], 1:7, type="n", xlab="AMCE", ylab="", axes=FALSE, xlim=c(-0.35,0.3), ylim=c(0.5,7), main= "New Leader")
for (i in 1:nrow(o4.r)){
	points(o4.r[i,1], nrow(o4.r)-(i-1), pch=16)
	lines(o4.r[i,c(2,3)], rep(nrow(o4.r)-(i-1),2))
	text(-0.3, nrow(o4.r)-(i-1), textLab[i], cex=0.8)
}
abline(v=0)
axis(1)

plot(o5.r[,1], 1:14, type="n", xlab="AMCE", ylab="", axes=FALSE, xlim=c(-0.35,0.3), ylim=c(0.5,14), main= "Military service")
for (i in 1:nrow(o5.r)){
	points(o5.r[i,1], nrow(o5.r)-(i-1), pch=16)
	lines(o5.r[i,c(2,3)], rep(nrow(o5.r)-(i-1),2))
	text(-0.3, nrow(o5.r)-(i-1), textLab[ifelse(i<8,i,i - 7)], cex=0.8)
}
abline(v=0)
axis(1)

plot(o6.r[,1], 1:7, type="n", xlab="AMCE", ylab="", axes=FALSE, xlim=c(-0.35,0.3), ylim=c(0.5,7), main= "Male Leader")
for (i in 1:nrow(o6.r)){
	points(o6.r[i,1], nrow(o6.r)-(i-1), pch=16)
	lines(o6.r[i,c(2,3)], rep(nrow(o6.r)-(i-1),2))
	text(-0.3, nrow(o6.r)-(i-1), textLab[i], cex=0.8)
}
abline(v=0)
axis(1)

plot(o7.r[,1], 1:7, type="n", xlab="AMCE", ylab="", axes=FALSE, xlim=c(-0.35,0.3), ylim=c(0.5,7), main= "Relationship with US")
for (i in 1:nrow(o7.r)){
	points(o7.r[i,1], nrow(o7.r)-(i-1), pch=16)
	lines(o7.r[i,c(2,3)], rep(nrow(o7.r)-(i-1),2))
	text(-0.3, nrow(o7.r)-(i-1), textLab[i], cex=0.8)
}
abline(v=0)
axis(1)

plot(o8.r[,1], 1:7, type="n", xlab="AMCE", ylab="", axes=FALSE, xlim=c(-0.35,0.3), ylim=c(0.5,7), main= "Previous opponent")
for (i in 1:nrow(o8.r)){
	points(o8.r[i,1], nrow(o8.r)-(i-1), pch=16)
	lines(o8.r[i,c(2,3)], rep(nrow(o8.r)-(i-1),2))
	text(-0.3, nrow(o8.r)-(i-1), textLab[i], cex=0.8)
}
abline(v=0)
axis(1)

plot(o9.r[,1], 1:7, type="n", xlab="AMCE", ylab="", axes=FALSE, xlim=c(-0.35,0.3), ylim=c(0.5,7), main= "Previous role")
for (i in 1:nrow(o9.r)){
	points(o9.r[i,1], nrow(o9.r)-(i-1), pch=16)
	lines(o9.r[i,c(2,3)], rep(nrow(o9.r)-(i-1),2))
	text(-0.3, nrow(o9.r)-(i-1), textLab[i], cex=0.8)
}
abline(v=0)
axis(1)

plot(o10.r[,1], 1:7, type="n", xlab="AMCE", ylab="", axes=FALSE, xlim=c(-0.35,0.3), ylim=c(0.5,7), main= "Previous action")
for (i in 1:nrow(o10.r)){
	points(o10.r[i,1], nrow(o10.r)-(i-1), pch=16)
	lines(o10.r[i,c(2,3)], rep(nrow(o10.r)-(i-1),2))
	text(-0.3, nrow(o10.r)-(i-1), textLab[i], cex=0.8)
}
abline(v=0)
axis(1)

plot(o11.r[,1], 1:7, type="n", xlab="AMCE", ylab="", axes=FALSE, xlim=c(-0.35,0.3), ylim=c(0.5,7), main= "Previous leader")
for (i in 1:nrow(o11.r)){
	points(o11.r[i,1], nrow(o11.r)-(i-1), pch=16)
	lines(o11.r[i,c(2,3)], rep(nrow(o11.r)-(i-1),2))
	text(-0.3, nrow(o11.r)-(i-1), textLab[i], cex=0.8)
}
abline(v=0)
axis(1)

plot(o12.r[,1], 1:14, type="n", xlab="AMCE", ylab="", axes=FALSE, xlim=c(-0.35,0.3), ylim=c(0.5,14), main= "Current behavior")
for (i in 1:nrow(o12.r)){
	points(o12.r[i,1], nrow(o12.r)-(i-1), pch=16)
	lines(o12.r[i,c(2,3)], rep(nrow(o12.r)-(i-1),2))
	text(-0.3, nrow(o12.r)-(i-1), textLab[ifelse(i<8,i,i - 7)], cex=0.8)
}
abline(v=0)
axis(1)

##### Main and weighted results in tabular form 

set.seed(43215)

#Main regression results
bootRegMod <- function(...){
	temp <- clusterBootS2(dat3)
	mod.temp <- lm(Chosen ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=temp)
	return(c(coef(mod.temp), sqrt(diag(vcov(mod.temp)))))
}

clusterExport(cl, list("dat3", "bootRegMod", "clusterBootS2"))

system.time(boot.regMod <- parSapply(cl, rep(1,1500), bootRegMod)) 

#Unweighted regression results
bootRegModW <- function(...){
	temp <- clusterBootS2(dat3)
	mod.temp <- lm(Chosen ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, weight=eWeight, data=temp)
	return(c(coef(mod.temp), sqrt(diag(vcov(mod.temp)))))
}

clusterExport(cl, list("bootRegModW"))

system.time(boot.regModW <- parSapply(cl, rep(1,1500), bootRegModW)) 

##### Table 3: Regression model specification of main results
## Tweak column labels in .tex
out1 <- cbind(round(apply(boot.regMod,1,mean)[1:17], digits=3), paste("(",round(apply(boot.regMod,1,mean), digits=3)[18:34],")", sep=""), paste("(",round(apply(boot.regMod,1,quantile,0.025), digits=3)[1:17], ", ", round(apply(boot.regMod,1,quantile,0.975), digits=3)[1:17], ")", sep=""))
out2 <- cbind(round(apply(boot.regModW,1,mean)[1:17], digits=3), paste("(",round(apply(boot.regModW,1,mean), digits=3)[18:34],")", sep=""), paste("(",round(apply(boot.regModW,1,quantile,0.025), digits=3)[1:17], ", ", round(apply(boot.regModW,1,quantile,0.975), digits=3)[1:17], ")", sep=""))
colnames(out1) <- colnames(out2) <- c("B", "SE", "95% CI")
rownames(out1) <- rownames(out2) <- c("(Intercept)", "High Capabilities", "High Stakes", "Democracy", "Mixed regime", "Adversary", "New Leader", "Long Military Service", "Some Military Service", "Male Leader", "Initiator", "Against adversary", "Different leader, backed down", "Same leader, backed down", "Same leader, stood firm", "Mobilized troops", "Public threat")
stargazer(cbind(out1, out2))

###### Results including the US
set.seed(43215)

bootConjoint <- function(...){
	temp <- clusterBootS2(dat2)
	mod.temp <- lm(Chosen ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("dat2", "bootConjoint", "clusterBootS2"))

system.time(boot.US <- parSapply(cl, rep(1,1500), bootConjoint)) 

#Extract quantities of interest
plot.US <- cbind(apply(boot.US, 1, mean), apply(boot.US, 1, quantile, c(0.025)), apply(boot.US, 1, quantile, c(0.05)), apply(boot.US, 1, quantile, c(0.95)), apply(boot.US,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.USb <- rbind(rep(0,5), plot.US[1,], rep(0,5), plot.US[2,], rep(0,5), plot.US[3:4,], rep(0,5), plot.US[5:6,], rep(0,5), plot.US[7,], rep(0,5), plot.US[c(9,8),], rep(0,5), plot.US[10,], rep(0,5), plot.US[11,], rep(0,5), plot.US[12,],rep(0,5), plot.US[c(15,14,13),], rep(0,5),plot.US[16:17,])

textLab <- c("Low Capabilities", "High Capabilities", "Low Stakes", "High Stakes", "Dictatorship", "Democracy", "Mixed", "Ally", "Adversary", "USA", "Established Leader", "New Leader", "No Military Service", "Some Military Service", "Long Military Service", "Female Leader", "Male Leader",  "Target", "Initiator", "Against ally", "Against adversary", "Different leader, stood firm", "Same leader, stood firm", "Same leader, backed down", "Different leader, backed down", "Nothing", "Mobilized troops", "Public threat")

#### Figure 4: Results including the US as a participant
# Clean up labels in post-production
dev.new(height=9, width=6.5)
par(mar=c(5.1,1.1,4.1,0), oma=c(0,4,0,0))
plot(plot.USb[,1], nrow(plot.USb):1, type="n", axes=FALSE, xlab="Average marginal component effect (AMCE)", ylab="",xlim=c(-0.3,0.2))

for (i in 1:nrow(plot.USb)){
	points(plot.USb[i,1], nrow(plot.USb)-(i-1), pch=16)
	lines(plot.USb[i,c(2,5)], rep(nrow(plot.USb)-(i-1),2))
	lines(plot.USb[i,c(3,4)], rep(nrow(plot.USb)-(i-1),2), lwd=2)		
	text(-0.3, nrow(plot.USb)-(i-1), textLab[i], cex=0.7)			
	}
abline(v=0)
axis(1)

###### Dropping unusual dyadic combinations
## Democratic dyads, allied US dyads ##
set.seed(43215)

#Find dyadic combinations to drop
dropList <- NULL
i <- unique(dat3$ID_numeric) 
for (j in i){
	for (r in 1:8){
		if (any(dat3$ID_numeric==j & dat3$Round==r)){
			locA <- which(dat3$ID_numeric==j & dat3$Round==r & dat3$Country=="A")
			locB <- which(dat3$ID_numeric==j & dat3$Round==r & dat3$Country=="B")
			#Any democratic dyads?
			if (dat3$RegimeType[locA]=="Democracy" & dat3$RegimeType[locB]=="Democracy"){
				dropList <- c(dropList, locA, locB)
			}
			#Any allied US dyads?
			if (dat3$actor[locA]=="Ally" & dat3$actor[locB]=="Ally"){
				dropList <- c(dropList, locA, locB)
			}			
		}
	}	
}

#How many observations are we dropping? 5788
length(dropList)

#Create new version of the dataset
dat.drop <- dat3[-dropList,]

#Re-estimating main model on pruned dataset
bootConjoint <- function(...){
	temp <- clusterBootS2(dat.drop)
	mod.temp <- lm(Chosen ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("dat.drop", "bootConjoint", "clusterBootS2"))

system.time(boot.drop <- parSapply(cl, rep(1,1500), bootConjoint)) 

#Extract quantities of interest
plot.matDrop <- cbind(apply(boot.drop, 1, mean), apply(boot.drop, 1, quantile, c(0.025)), apply(boot.drop, 1, quantile, c(0.05)), apply(boot.drop, 1, quantile, c(0.95)), apply(boot.drop,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.matDropb <- rbind(rep(0,5), plot.matDrop[1,], rep(0,5), plot.matDrop[2,], rep(0,5), plot.matDrop[3:4,], rep(0,5), plot.matDrop[5,], rep(0,5), plot.matDrop[6,], rep(0,5), plot.matDrop[c(8,7),], rep(0,5), plot.matDrop[9,], rep(0,5), plot.matDrop[10,], rep(0,5), plot.matDrop[11,],rep(0,5), plot.matDrop[c(14,13,12),], rep(0,5),plot.matDrop[15:16,])

###If you haven't run lines 35-39 and 48-57 from BJPS Replication 2.R, run the code below:

set.seed(43215)

bootConjoint <- function(...){
	temp <- clusterBootS2(dat3)
	mod.temp <- lm(Chosen ~  Capabilities + Stakes + RegimeType + actor + NewLeader + milService + maleLeader + prevRole + prevOpponent + prevActLdr + curBehavior, data=temp)
	return(coef(mod.temp))
}

clusterExport(cl, list("dat3", "bootConjoint", "clusterBootS2"))

#Main analysis
system.time(boot.full2 <- parSapply(cl, rep(1,1500), bootConjoint)) #Main results

#Extract quantities of interest
plot.mat2 <- cbind(apply(boot.full2, 1, mean), apply(boot.full2, 1, quantile, c(0.025)), apply(boot.full2, 1, quantile, c(0.05)), apply(boot.full2, 1, quantile, c(0.95)), apply(boot.full2,1,quantile, c(0.975)))[-1,]

#Now add baseline categories
plot.mat2b <- rbind(rep(0,5), plot.mat2[1,], rep(0,5), plot.mat2[2,], rep(0,5), plot.mat2[3:4,], rep(0,5), plot.mat2[5,], rep(0,5), plot.mat2[6,], rep(0,5), plot.mat2[c(8,7),], rep(0,5), plot.mat2[9,], rep(0,5), plot.mat2[10,], rep(0,5), plot.mat2[11,],rep(0,5), plot.mat2[c(14,13,12),], rep(0,5),plot.mat2[15:16,])

######

textLab <- c("Low Capabilities", "High Capabilities", "Low Stakes", "High Stakes", "Dictatorship", "Democracy", "Mixed", "Ally", "Adversary", "Established Leader", "New Leader", "No Military Service", "Some Military Service", "Long Military Service", "Female Leader", "Male Leader",  "Target", "Initiator", "Against ally", "Against adversary", "Different leader, stood firm", "Same leader, stood firm", "Same leader, backed down", "Different leader, backed down", "Nothing", "Mobilized troops", "Public threat")

##### Figure 5: Robustness check: dropping unusual combinations
## Clean up labels in post-production
dev.new(height=9, width=6.5)
par(mar=c(5.1,1.1,4.1,0), oma=c(0,4,0,0))
plot(plot.matDropb[,1], nrow(plot.matDropb):1, type="n", axes=FALSE, xlab="Average marginal component effect (AMCE)", ylab="",xlim=c(-0.3,0.2))

for (i in 1:nrow(plot.matDropb)){
	points(plot.mat2b[i,1], nrow(plot.mat2b)-(i-1), pch=16)
	points(plot.matDropb[i,1], nrow(plot.matDropb)-(i-1)-0.25, pch=15, col="grey")
	lines(plot.mat2b[i,c(2,5)], rep(nrow(plot.mat2b)-(i-1),2))
	lines(plot.matDropb[i,c(2,5)], rep(nrow(plot.matDropb)-(i-1)-0.25,2), col="grey")
	lines(plot.mat2b[i,c(3,4)], rep(nrow(plot.mat2b)-(i-1),2), lwd=2)		
	lines(plot.matDropb[i,c(3,4)], rep(nrow(plot.matDropb)-(i-1)-0.25,2), lwd=2, col="grey")		
	text(-0.3, nrow(plot.matDropb)-(i-1), textLab[i], cex=0.7)			
	}
abline(v=0)
axis(1)
